From GPA to Shape Data

• In this lecture we:
  • Obtain variables we can use for analyses
  • Review the relative position of individuals in shape space
  • Visualize shape variation

Shape Spaces from GPA

• After GPA, each landmark configuration is a point in shape space
• The ‘center’ of this universe is the consensus (global mean shape)
• Each shape occupies a unique point in shape space
• The metric of shape space is the Procrustes distance \(\small{D}_{Proc}\)

\[\tiny{D}_{Proc}=\sqrt{\sum_{i,j}^{k,p}\left(\mathbf{Y}_{1.ij}-\mathbf{Y}_{2.ij}\right)^2}\]

Shape Spaces from GPA

• After GPA, each landmark configuration is a point in shape space
• The ‘center’ of this universe is the consensus (global mean shape)
• Each shape occupies a unique point in shape space
• The metric of shape space is the Procrustes distance \(\small{D}_{Proc}\)

\[\tiny{D}_{Proc}=\sqrt{\sum_{i,j}^{k,p}\left(\mathbf{Y}_{1.ij}-\mathbf{Y}_{2.ij}\right)^2}\]

• Shape space is curved!

Example with 2,000 random (uniform) triangles. Note overabundance of shapes near ‘north pole’.

Consequences of GPA

• With GPA, specimen location, scale and rotation are standardized, resulting in redundant dimensions

• \(\small{2p-4}\) for 2D landmarks or \(\small{3p-7}\) for 3D landmarks

• Statistical Challenges:
  • Covariance matrices are singular, so statistical hypothesis testing via parametric approaches cannot be completed
  • Shape space is curved, and most standard statistics assume linear spaces

Consequences of GPA

• With GPA, specimen location, scale and rotation are standardized, resulting in redundant dimensions

• \(\small{2p-4}\) for 2D landmarks or \(\small{3p-7}\) for 3D landmarks

• Statistical Challenges:
  • Covariance matrices are singular, so statistical hypothesis testing via parametric approaches cannot be completed
  • Shape space is curved, and most standard statistics assume linear spaces
• Solutions:
  • Data projection to a tangent space. This provides both a linear approximation to shape space and allows for redudant dimensions to be removed
  • Use RRPP (residual randomization in permutation procedures) for evaluating statistical hypotheses (subsequent lectures)

Kendall’s Tangent Space Coordinates

• Orthogonal projection of shapes to linear tangent space

\[\small\mathbf{X}'=\mathbf{X\left(I_{kp}-X_c^T(X_cX_c^T)^{-1}X_c\right)}\]

• Specimen scores from projection are Kendall’s tangent space coordinates (Dryden and Mardia, 1993; Kent 1994)
• Perform PCA and retain only 2p-4 (3D: 3p-7) shape variables with variation if desired (NOTE: permutation via RRPP is now standard practice, so this step is not strictly necessary!)

Kendall’s Tangent Space Coordinates

• Random uniform shapes ARE uniformally distributed in Tangent Space!

Same set of 2,000 randomly generated triangles

Shape Differences as Deformations

• With Tangent Space we can view patterns of shape variation; what about visualzing individual shape differences?
• D’Arcy Thompson, in “On Growth and Form” (1917) visualized shape differences using transformation grids
  
• Simple mathematical models to represent complex morphological changes

• Today we use the Thin-Plate Spline to generate such visualizations

Shape Deformations: the Thin-Plate Spline

• Based on D’Arcy Thompson’s, Bookstein (1989) proposed using the thin-plate spline for providing a mathematically accurate representation of shape change (math was borrowed from engineering methods from 1980s)
• Deformation of a reference to a target specimen
• The coefficients of the fitted TPS model describe shape changes

• The TPS is an alternative way of projecting aligned shape coordinates to tangent space (different math, but same result: shown later)

The Thin-Plate Spline

• TPS is a smooth interpolation function (i.e., a doubly-differentiated equation)
• It models the differences in landmark locations between the reference and target specimen for each coordinate dimension separately (x,y)

The Thin-Plate Spline

• TPS is a smooth interpolation function (i.e., a doubly-differentiated equation)
• It models the differences in landmark locations between the reference and target specimen for each coordinate dimension separately (x,y)

General TPS Model for \(\small{1}\) landmark: \[\tiny\begin{bmatrix} x^1\\y^1\end{bmatrix}=A\begin{bmatrix}1\\x\\y\end{bmatrix}+\sum_{i=1}^pw_iU(r_i)\]

where \(\tiny{A}\begin{bmatrix}1\\x\\y\end{bmatrix}\) represents the affine (uniform) component and \(\small\sum_{i=1}^pw_iU(r_i)\) represents to non-affine (non-uniform) component

TPS Computations

• 1: From the reference (consensus) estimate: \(\tiny\mathbf{L}= \left[ \begin{array}{c|c} \mathbf{P} & \mathbf{Q}\\ \hline \mathbf{Q}^T & \mathbf{0} \\ \end{array} \right]\)

where: \(\tiny\mathbf{Q}= \begin{bmatrix} 1 & x_1 & y_1 \\ 1 & x_2 & y_2 \\ \vdots & \vdots & \vdots \\ 1 & x_p & y_p \end{bmatrix}\) & \(\tiny\mathbf{0}= \begin{bmatrix} 0 & 0 & 0 \\ 0 & 0 & 0 \\ 0 & 0 & 0 \end{bmatrix}\)

\(\tiny\mathbf{P}=\tiny\begin{bmatrix} 0 & U(r_{12}) & U(r_{12}) & \dots & U(r_{1p}) \\ U(r_{21}) & 0 & U(r_{23}) & \dots & U(r_{2p}) \\ U(r_{31}) & U(r_{32}) & 0 & \dots & U(r_{3p}) \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ U(r_{p1}) & U(r_{p2}) & U(r_{p3}) & \dots & 0 \end{bmatrix}\)

TPS Computations

• 1: From the reference (consensus) estimate: \(\tiny\mathbf{L}= \left[ \begin{array}{c|c} \mathbf{P} & \mathbf{Q}\\ \hline \mathbf{Q}^T & \mathbf{0} \\ \end{array} \right]\)

where: \(\tiny\mathbf{Q}= \begin{bmatrix} 1 & x_1 & y_1 \\ 1 & x_2 & y_2 \\ \vdots & \vdots & \vdots \\ 1 & x_p & y_p \end{bmatrix}\) & \(\tiny\mathbf{0}= \begin{bmatrix} 0 & 0 & 0 \\ 0 & 0 & 0 \\ 0 & 0 & 0 \end{bmatrix}\)

\(\tiny\mathbf{P}=\tiny\begin{bmatrix} 0 & U(r_{12}) & U(r_{12}) & \dots & U(r_{1p}) \\ U(r_{21}) & 0 & U(r_{23}) & \dots & U(r_{2p}) \\ U(r_{31}) & U(r_{32}) & 0 & \dots & U(r_{3p}) \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ U(r_{p1}) & U(r_{p2}) & U(r_{p3}) & \dots & 0 \end{bmatrix}\)

Here, \(\small\mathbf{Q}\) contains the coordinates of the reference and \(\small\mathbf{P}\) is found as: \(\small{U(r_{ij})}=r_{ij}^2ln(r_{ij})\) with \(\small{r_{ij}}\) as the distance between the \(\small{i^{th}}\) and \(\small{j^{th}}\) landmarks

NOTES: \(\small{0}\) is \(\small{4 \times 4}\) for 3D data. Also, for 3D data: \(\small{U(r_{ij})}=|{r_{ij}}|\)

TPS Computations (Cont.)

• 2: Obtain \(\small\mathbf{L}^{-1}\). Upper-left \(\small{p \times p}\) block is the Bending Energy Matrix (\(\small\mathbf{L}_p^{-1}\))

TPS Computations (Cont.)

• 2: Obtain \(\small\mathbf{L}^{-1}\). Upper-left \(\small{p \times p}\) block is the Bending Energy Matrix (\(\small\mathbf{L}_p^{-1}\))
• 3: Perform eigenanalysis of \(\small\mathbf{L}_p^{-1}\) to obtain Principal Warps which are the \(\small{p-3}\) eigenvectors containing variation (\(\small{p-4}\) for 3D)

\[\small\mathbf{L}_p^{-1}=\mathbf{E\Lambda{E}^T}\]

TPS Computations (Cont.)

• 2: Obtain \(\small\mathbf{L}^{-1}\). Upper-left \(\small{p \times p}\) block is the Bending Energy Matrix (\(\small\mathbf{L}_p^{-1}\))
• 3: Perform eigenanalysis of \(\small\mathbf{L}_p^{-1}\) to obtain Principal Warps which are the \(\small{p-3}\) eigenvectors containing variation (\(\small{p-4}\) for 3D)

\[\small\mathbf{L}_p^{-1}=\mathbf{E\Lambda{E}^T}\]

• 4: Estimate non-uniform shape variables (Partial Warp Scores) via projection onto principal warps \(\small\mathbf{E}\):

\[\small\mathbf{W}=\mathbf{V(I_2\otimes{E})=V}\begin{bmatrix}\mathbf{E} & \mathbf{0} \\ \mathbf{0} & \mathbf{E} \end{bmatrix})\]

where \(\small\mathbf{V}=\mathbf{[V_x|V_y]}\) with \(\small\mathbf{V_x}=\mathbf{(X-\overline{X}_{ref})}\) & \(\small\mathbf{V_y}=\mathbf{(Y-\overline{Y}_{ref})}\)

TPS Computations (Cont.)

• 2: Obtain \(\small\mathbf{L}^{-1}\). Upper-left \(\small{p \times p}\) block is the Bending Energy Matrix (\(\small\mathbf{L}_p^{-1}\))
• 3: Perform eigenanalysis of \(\small\mathbf{L}_p^{-1}\) to obtain Principal Warps which are the \(\small{p-3}\) eigenvectors containing variation (\(\small{p-4}\) for 3D)

\[\small\mathbf{L}_p^{-1}=\mathbf{E\Lambda{E}^T}\]

• 4: Estimate non-uniform shape variables (Partial Warp Scores) via projection onto principal warps \(\small\mathbf{E}\):

\[\small\mathbf{W}=\mathbf{V(I_2\otimes{E})=V}\begin{bmatrix}\mathbf{E} & \mathbf{0} \\ \mathbf{0} & \mathbf{E} \end{bmatrix})\]

where \(\small\mathbf{V}=\mathbf{[V_x|V_y]}\) with \(\small\mathbf{V_x}=\mathbf{(X-\overline{X}_{ref})}\) & \(\small\mathbf{V_y}=\mathbf{(Y-\overline{Y}_{ref})}\)

Note, one can weight partial warp scores inversely by spatial scale as: \(\small\mathbf{W}=\mathbf{V(I_2\otimes{E\Lambda^{-\alpha/2}})}\)

• There are \(\small{2p-6}\) partial warp scores (\(\small{3p-12}\) for 3D data)

TPS Model: Uniform (Affine) Shape Estimation

• Affine shape describes shape changes where parallel lines remain parallel
• Affine shape found as the complement to non-affine
• Total shape (from GPA) is the sum of affine and non-affine shape:

\[\small\mathbf{S=U\oplus{B}}\]

TPS Model: Uniform (Affine) Shape Estimation

• Affine shape describes shape changes where parallel lines remain parallel
• Affine shape found as the complement to non-affine
• Total shape (from GPA) is the sum of affine and non-affine shape:

\[\small\mathbf{S=U\oplus{B}}\]

• Computations
  • Construct space orthogonal to non-affine shape component: \(\small\mathbf{N=I_p-E(E^TE)^{-1}E^T}\)
  • Project specimens: \(\small\mathbf{V(N\oplus{I_k})}\)
  • Extract affine (uniform) dimensions through singular-value decomposition (SVD): \(\small\mathbf{LSR^T=V(N\oplus{I_k})}\). There are \(\small{2}\) uniform dimensions (\(\small{5}\) for 3D)

Rohlf and Bookstein. (2003). Syst. Biol.

Total Shape Variation

• Shape variation includes BOTH uniform and non-uniform components

• Together, the Partial Warps + Uniform component scores comprise the total set of shape variables: \(\small\mathbf{W}=\mathbf{V(I_2\otimes{E})|U}\)

• There are \(\small{2p-4}\) shape variables for 2D data and \(\small{3p-7}\) for 3D data

TPS Model: Example

• Model shape differences from a square to a kite (only non-uniform shape differences present)

\(\tiny\mathbf{X}_{sqr}=\begin{bmatrix} 0 & 0.5\\ -0.5 & 0\\ 0 & -0.5\\ 0.5 & 0 \end{bmatrix}\) and \(\tiny\mathbf{X}_{kite}=\begin{bmatrix} 0 & 0.375\\ -0.5 & 0.125\\ 0 & -0.625\\ 0.5 & 0.125 \end{bmatrix}\)

TPS Model: Example (Cont.)

• Obtain: \(\tiny\mathbf{L}= \left[ \begin{array}{c|c} \mathbf{P} & \mathbf{Q}\\ \hline \mathbf{Q}^T & \mathbf{0} \\ \end{array} \right]=\left[ \begin{array}{cccc|ccc} 0 & -0.173 & 0 & -0.173 & 1 & 0 & 0.5\\ -0.173 & 0 & -0.173 & 0 & 1 & -0.5 & 0 \\ 0 & -0.173 & 0 & -0.173 & 1 & 0 & -0.5\\ -0.173 & 0 & -0.173 & 0 & 1 & -0.5 & 0 \\ \hline 1 & 1 & 1 & 1 & 0 & 0 & 0 \\ 0 & -0.5 & 0 & 0.5 & 0 & 0 & 0 \\ 0.5 & 0 & -0.5 & 0 & 0 & 0 & 0 \\ \end{array} \right]\)

• Then calculate: \(\tiny\mathbf{L^{-1}}=\left[ \begin{array}{cccc|ccc} 0.7271 & -0.7271 & 0.7271 & -0.7271 & 0.25 & 0 & 1\\ -0.7271 & 0.7271 & -0.7271 & 0.7271 & 0.25 & -1 & 0 \\ 0.7271 & -0.7271 & 0.7271 & -0.7271 & 0.25 & 0 & -1\\ -0.7271 & 0.7271 & -0.7271 & 0.7271 & 0.25 & 1 & 0 \\ \hline 0.25 & 0.25 & 0.25 & 0.25 & 0.087 & 0 & 0 \\ 0 & -1 & 0 & 1 & 0 & 0 & 0 \\ 1 & 0 & -1 & 0 & 0 & 0 & 0 \\ \end{array} \right]\)

• \(\tiny\mathbf{L_p^{-1}}=\left[ \begin{array} 0.7271 & -0.7271 & 0.7271 & -0.7271 \\ -0.7271 & 0.7271 & -0.7271 & 0.7271 \\ 0.7271 & -0.7271 & 0.7271 & -0.7271 \\ -0.7271 & 0.7271 & -0.7271 & 0.7271 \\ \end{array} \right]\)

TPS Model: Example (Cont.)

• Solve: \(\small\mathbf{L}_p^{-1}=\mathbf{E\Lambda{E}^T}\)

\(\tiny\mathbf{L_p^{-1}}=\left[ \begin{array} 0.7271 & -0.7271 & 0.7271 & -0.7271 \\ -0.7271 & 0.7271 & -0.7271 & 0.7271 \\ 0.7271 & -0.7271 & 0.7271 & -0.7271 \\ -0.7271 & 0.7271 & -0.7271 & 0.7271 \\ \end{array} \right]=\)

• \(\tiny\mathbf{E}=\left[ \begin{array} \mathbf{0.5} & 0.866 & 0 & 0 \\ \mathbf{-0.5} & 0.289 & 0 & -0.0816 \\ \mathbf{0.5} & -0.289 & 0.707 & -0.408 \\ \mathbf{-0.5} & 0.289 & -0.707 & 0.408 \\ \end{array} \right]\) & \(\tiny\mathbf{\Lambda}=\left[ \begin{array} \mathbf{2.885} & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 \\ \end{array} \right]\)

• Obtain shape variables: \(\small\mathbf{W=VE}\) using: \(\tiny\mathbf{V_{kite}}=\begin{bmatrix} 0 & 0 & 0 & 0 & -1.25 & -1.25 & -1.25 & -1.25 \end{bmatrix}\)

• \(\tiny\mathbf{W=VE}=\mathbf{V_{kite}}\begin{bmatrix} 0.5 & 0 \\ -0.5 & 0 \\ 0.5 & 0 \\ -0.5 & 0 \\ 0 & 0.5 \\ 0 & -0.5 \\ 0 & 0.5 \\ 0 & -0.5 \\ \end{bmatrix}=\begin{bmatrix} 0 \\ -0.25 \end{bmatrix}\)

• NOTE: Uniform shape variables, \(\small\mathbf{U_1}=0\) & \(\small\mathbf{U_2} = 0\) for this example

Which Shape Variables to Use?

• We now have 3 shape variable options:
  • GPA-aligned coordinates
  • Kendall’s Tangent Space coordinates (GPA + orthogonal projection)
  • TPS variables (partial warp scores + uniform shape)
• Which ones should we use?

Which Shape Variables to Use?

• We now have 3 shape variable options:
  • GPA-aligned coordinates
  • Kendall’s Tangent Space coordinates (GPA + orthogonal projection)
  • TPS variables (partial warp scores + uniform shape)
• Which ones should we use?

• Plot of shape distances for Podarcis data obtained from GPA-shape space and Kendall’s Tangent space. Correlation of 1.0 in this case (usually very, very high)

• AND correlation between PWS+U vs. Tangent space is ALWAYS 1.0 (see Rohlf 1999)

Shape Variables: Conclusions

• Proc. Tangent coordinates & TPS shape variables yield identical results
  • They are simply rotations of one another
• GPA-aligned coordinates (unprojected) are close, but not exact
  • Difference due to curvature of shape space

• Use Kendall´s Tangent Space Coordinates to represent shape
• These are found from GPA+Orthogonal Projection

• NOTE: Most current papers use tangent space coordinates, aka Procrustes residuals combined with resampling tests (to accomodate singular dimensions)
• TPS then, is frequently used to visualize shape differences

Exploring Shape Variation

• GM data are points in high-dimensional shape space
  • (e.g. for \(\small{p=20}\) & \(\small{k=3}\), shape space comprises 53 dimensions)
• Our human brains are limited to perception in up to 3 dimensions
• We need tools for visualizing and exploring shape variation in high-dimensional data spaces

• Principal Components Analysis (PCA) is one such tool

Principal Components Analysis (PCA)

• PCA combines orthogonal rotation and projection to arrive at a summary plot of the dataspace
• The objective: to summarize most of the variation in as few dimensions as possible
• It consists of a rigid rotation of the data based on directions of variation, followed by projection to those new summary axes
• First one finds the set of axes the describe progressively less variation in the data
• Next, the data are rotated so that the main axis of variation (PC1) is horizontal
• Subsequent axes are orthogonal to PC1

PCA: Conceptual Visualization

• This is what PCA does:

PCA: Conceptual Visualization

• This is what PCA does:

• PCA has performed a rigid rotation of the original data, so that variation is aligned with the new (PC) axes

PCA: Standard Computations

Principal component analysis (PCA) is two things: (1) a singular-value decomposition (SVD) of a symmetric matrix and (2) projection of mean-centered or standardized data onto eigenvectors from SVD.

Using mean-centered data: \(\small\mathbf{Y}_c\) we calculate:

\[\small\hat{\mathbf{\Sigma}}=\frac{\mathbf{Y}^{T}_c\mathbf{Y}_c}{n-1}\]

We then decompose the covariance matrix via eigen-analysis (SVD):

\[\small\hat{\mathbf{\Sigma}}=\mathbf{U} \mathbf{\Lambda} \mathbf{U}^T\]

This step identifies a set of orthogonal axes describing directions of variation. These are the columns of \(\small\mathbf{U}\). Next we project the data onto these axes as:

\[\small\mathbf{P} = \mathbf{Y}_c\mathbf{U}\]

Here, \(\small\mathbf{P}\) represent the set of projection scores (principal component scores) on each of the PC axes. These are used in the subsequent plot.

Finally, the percent variation explained by any principal component (column of \(\mathbf{U}\)) is \[\frac{\lambda_i}{tr\mathbf{\Lambda}}\]

PCA: Alternative Computations

Note that one can alternatively perform singular-value decomposition (SVD) directly on \(\small\mathbf{Y}_c\):

\[\small{\mathbf{Y}_c}=\mathbf{V} \mathbf{D} \mathbf{U}^T\]

Here, \(\small\mathbf{D^2}\) expresses the percent-variation explained by each PC-axis, which are found as the right-singular vectors in \(\small\mathbf{U}\).

PC scores are found as:

\[\small\mathbf{P} = \mathbf{VD}\]

NOTE: eigenanalysis of the inner product \(\small\mathbf{Y}^{T}_c\mathbf{Y}_c\) yields \(\small\mathbf{U}\), while eigenanalysis of the outer product \(\small\mathbf{Y}_c\mathbf{Y}^{T}_c\) yields \(\small\mathbf{V}\)

PCA of Shape Data: Relative Warps

• SVD of shape variation

\[\small{\mathbf{W^*}}=\mathbf{V} \mathbf{D} \mathbf{U}^T\]

• Here, \(\small\mathbf{VD}\) are the PC scores (relative warp scores) of the shape data

PCA of Shape Data: Relative Warps

• SVD of shape variation

\[\small{\mathbf{W^{*}}}=\mathbf{V} \mathbf{D} \mathbf{U}^T\]

• Here, \(\small\mathbf{VD}\) are the PC scores (relative warp scores) of the shape data

• NOTE: if W found as: \(\small\mathbf{W}=\mathbf{V(I_2\otimes{E\Lambda^{-\alpha/2}})}\) PCA=RWA is weighted by spatial scale of deformations

• In this case
  • \(\small\alpha = 0\): unweighted analysis
  • \(\small\alpha < 0\): emphasizes small-scale variation
  • \(\small\alpha > 0\): emphasizes large-scale variation

• Careful! \(\small\alpha\neq{0}\) distorts data space, so it changes the relationships between observations!

RWA: Example

PCA of Podarcis data

RWA: Example (Cont.)

• PCA using different \(\small\alpha\) values (\(\small{1}\) & \(\small{-1}\))

RWA: Thoughts

• Ordination through PCA very useful approach for exploring shape space

• \(\small\alpha\neq{0}\) in RWA allows emphasizing small- or large-scale deformations

• However, \(\small\alpha\neq{0}\) changes shape space, so we need to choose carefully

• In practice, there is generally no conceptual reason of biological relevance for emphasizing deformations at different spatial scales

• \(\small\alpha{=0}\) is preferred: no a priori, arbitrary emphasis on specific variables, and it allows the use of total shape space (i.e. PW scores + U)

Visualizing Shape Differences

• TPS allows us to visualize what exactly changes in shape

• As an interpolation function, TPS allows us to precisely map all points of one shape to those of another shape, and visualize how they differ

• One may wish to visualize:
  • Differences between groups
  • Extremes of PC axes
  • How shape changes with size (allometry)
  • Shape evolution (differences between ancestor and descendent on a phylogeny)

Visualizing Shape Differences

• TPS allows us to visualize what exactly changes in shape

• As an interpolation function, TPS allows us to precisely map all points of one shape to those of another shape, and visualize how they differ

• One may wish to visualize:
  • Differences between groups
  • Extremes of PC axes
  • How shape changes with size (allometry)
  • Shape evolution (differences between ancestor and descendent on a phylogeny)

• One can also magnify shape differences to aid biological interpretation

*TPS of Ref->Target is a vector in shape space, which can be extended in length

Visualizing Shape Differences: Example

• TPS of male vs. female Podarcis (5X magnification)

• Major difference in posterior region of skull

Shape Variables: Conclusions

• Produce shape variables through projection to tangent space:
  • Kendall´s tangent space coordinates (Procrustes residuals [GPA+Projection])
• Explore shape variation through ordination in tangent space:
  • PCA of shape data (aka Relative Warps)
• Visualize shape variation to make biological interpretations:
  • Use TPS to produce deformation grids